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ABSTRACT 

This work was motivated by the need that may arise 
during the creation of finite element programs, for 
higher order elements. The problem of selecting or 
creating shape functions that satisfy the required need 
of the problem may be one major problem that stand 
in the way of the element creation. In this work we 
present an attempt that points in the direction of 
creating the element matrices for higher order 
elements using simple, Lagrange, and modified 
Lagrange polynomials. The results obtained for the 
test cases indicate the possibilities and limitations on 
those attempts. It is concluded that the modifies 
Lagrange polynomials have a great potential for use in 
elements with high number of nodes and/or high 
number of DOL per node. Nevertheless, they impose a 
high computational cost on the element matrices 
generation. 

Keywords: higher order elements, finite element 
analysis 

1. INTRODUCTION 

In the past decades, the finite element method has 
grown into a well developed numerical technique that 
can be used to provide accurate and relatively quick 
solution to boundary value problems that describe 
several physical problems. The procedure of the finite 
element model starts by dividing the domain into 
elements that are connected through nodes. Each node 
has degrees of freedom (DOE) that, usually, present 
the value and derivatives of the physical function of 
interest. Then the elements are defined by the 
interpolation functions that describe the change of the 
function within their boundaries. Lollowing that, the 
whole domain description is created by assembling 
the elements and applying the boundary conditions. 
Linally, the problem is solved to obtain the values of 
the DOL of all the nodes. 1, 2 


The basic assumption is that as the number of 
elements increase in the domain, the solution 
approaches the exact one. The power of the finite 
element method lies in that it may use simple 
functions, usually polynomials, to describe the change 
of the target function inside an element, and then use 
that to describe the change in the whole domain by 
assembling the different elements that cover the 
domain. 

Another approach to reach higher accuracy is to 
increase the order of each of the elements used in the 
model which reduces the need for higher number of 
elements (p-version finite element). Higher order 
elements may be generated by using more nodes in 
the element, using more complex functions, or by 
hierarchical elements. 3, 4, 5, 6 

Most of the work performed for the solution of the 
finite element problem is algorithmic, that enabled the 
excessive use computers and the production of 
computer packages that can manipulate several 
physical problems with complex domain geometries. 
Nevertheless, in the core of all finite element codes, 
lies the element matrices. The element matrices are 
generated through the knowledge of the mathematical 
model that describes the physical problem and the 
selection of the interpolation functions, shape 
functions. Equipped with those two pieces of 
information, you may generate the element matrices 
for any physical problem. However, the generation of 
a general element is still far from algorithmic. 7 

A special problem faces the generation of the finite 
element model is the selection of the interpolation 
function that satisfies the prescribed degrees of 
freedom at the nodes 5, 6, 8. Many of the problems 
solved involved the selection of functions that are 
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published in the literature most of which are the so 
called Lagrange Polynomials. The Lagrange 
polynomials possess a preferred characteristic that 
they have the value of one at a specific point, node, in 
the domain and the value of zero at all other nodes. 
When the DOF involve derivatives of the function, 
the Lagrange polynomials should be selected to 
satisfy the condition that they have the slope of one at 
the associated node and zero at the others as well. 
However, it is quite difficult to find the Lagrange 
polynomials that may be used in a general problem 
that is why they may be generated by using classical 
interpolation techniques. 

The classical interpolation techniques require the 
evaluation of the polynomial coefficients by forcing 
the polynomial to satisfy the boundary conditions, 
thus, requiring the solution of a set algebraic equation, 
hence, the inversion of a matrix. The matrix inversion 
process is a straight forward numerical technique; 
however, as the number of variables increase, the 
matrix eventually becomes singular because of the 
round-off errors. 9, 10 

Another challenge that faces the creators of finite 
element models is the performance of integration. The 
element matrices are generated by integrating the 
shape functions or their derivatives over the element. 
In many simple cases, the integration is readily 
available by hand or using symbolic manipulators. 
However, in more practical cases, numerical 
integration is required which, in turn, introduces 
errors to the resulting element models. 

In the past years, several research articles were 
published, mostly by mathematics oriented 
researchers, about the topic of generating higher order 
polynomials, shape functions, for applications in the 
finite element models 5, 6, 8, 11, 12. Such research, 
among many others not cited here, accomplished a 
great task in the direction of generalizing the 
automation of finite element model generation. 
Following that, an excellent attempt for the 
automation was presented in 7 using the symbolic 
manipulator Mathematica® and an accompanying 
package AceGen®. 

None of the research reviewed by the author presented 
a straight-forward methodology that may be used by 
the engineers who need to create finite element 
models on numerical manipulators, thus, in the 
following work, we will be attempting to describe a 


generalize procedure for the generation of the finite 
element matrices with an eye on elements with large 
numbers of nodes and degrees of freedom. We will 
start by describing a classical method of generating 
the matrices and creating the shape functions, then we 
will move towards describing attempts to generate a 
general method for performing exact integration for 
generating element matrices. Following that, we will 
describe the attempts using Lagrange polynomials to 
avoid matrix inversion, finally, modifies Lagrange 
polynomials will be created to enable to the creation 
of different finite element problems. 

The main purpose of this work is providing a 
methodology that may by applied using numerical 
manipulators, especially open source, to generate the 
element matrices that may be used by researchers 
without the need for commercial packages or 
commercial symbolic manipulators. 

2. CLASSICAL APPROACH 

2.1. Interpolation Polynomials - Shape Functions 

The finite element model starts by assuming the 
solution of the unknown function in terms of an 
interpolation polynomial that satisfies the values of 
the degrees of freedom at the nodes of the element. In 
most text, these functions are denoted the letter N 
such that: 

f(x)= f,N,(x)+f 2 N 2 (x)+... 

Where / are the values of the function at node i and 
Ni(x) are polynomials that have the value of one at the 
corresponding node and zero at every other node. The 
above relation may be written as: 

f(x)=[N,(x) N 2 (x) ...\lfl\=[N(x)\{5} 


The first problem that faces the researcher who wants 
to create a new finite element model is to find and 
select the polynomials that may be used for the 
problem. However, that is not a great problem for low 
order elements since the polynomials are readily 
available in many texts. Nevertheless, coding them 
may involve mistyping which will require debugging 
of the code. On the other hand, it becomes a bit 
inconvenient to change the number of nodes per 
element to test the effect of higher, or lower, order 
elements. Thus we may use another way of defining 
the polynomials that can reduce both problems, that is 
by using regular simple polynomials such that: 
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f(x)=a 0 +aix+a 2 x z ...= ['\ x x 2 ...J 


3 | 


r[H(x)\{a} 


Where a, are constants, generalized coordinates that 
should be determined for the function to satisfy the 
given values. To do that we will need to set a number 
of equations and solve them to find those values, thus: 
f(x i)= f,= [H(x,)\{a] 
f(x 2 )= f 2 = [H(x 2 )\{a } 
f(x 3 h f 3 = [H(x 3 )\{a] 


Which may be written as: 

mm 



[H(x 2 )\ 

[H(x 3 )\ 


[T]{a}={5} 

Which gives: 

{a}=[r 1 ]{5} 
f(xh [H{x)\[T']{5} 

It can be readily proven that: 

LA/(*)j= [H(x)\[T'] 

Thus using this procedure we were able to obtain the 
polynomials for any number of degrees of freedom 
without having to look them up in literature or 
enduring the problems of mistyping them into the 
code. However, nothing comes for free; using this 
method will involve matrix inversion which will 
become ill-conditioned as the number of degrees of 
freedom increase. 

If the degrees of freedom of the element involve the 
derivatives of the function, we may adjust the above 
procedure. Given the value of the function/, and the 
fist derivative/’,, for example, we may write: 


f(x)=[H\x)\{a}=[0 1 2x ...J, 


Then the [T] matrix may be constructed such that: 

■L H(x,)V 


[Th 


[H'(x,)\ 

[H(x 2 )\ 

[H'(x 2 )\ 


In this case, the number of degrees of freedom per 
element will be equal to twice the number of nodes 


and the function N,(x) will appear in pairs. The first 
function of the pair will have a value of one at its 
corresponding node while the slope at that node is 
zero, and has a value and slope of zero at every other 
node. Meanwhile, the second function of the pair will 
have a slope of one and a value of zero at the 
corresponding node while its value and slope are 
equal to zero at every other node. 

Note that, in both cases above, the order of the 
polynomial terms in the H(x) matrix does not affect 
the resulting shape functions because they are going 
to be rearranged using the T matrix. 

2.2. Element Matrices 

Each finite element problem will generate a set of 
matrices for each element. These matrices are derived 
from the physical model that is usually presented in 
the form of a differential equation. The matrices are 
constructed by integrating some derivative of the 
function that usually appears in the form: 
r 2 

k = \ Q{x) D m {[N{x)\) T D m {[N{x)\)dx 

Where Q(x) is some function that describes the 
physical properties of the problem, D m (.) is a 
differential operator that differentiates the function 
f(x) m times. For example, in mechanics of material, 
the finite element matrix for a bar in static problem 
may be obtained using: 

E(x)A(x){N'(x)}[N'(*)ldx 

Where E(x) and A(x) are the modulus of elasticity and 
the cross-section area respectively. For simple 
problems, the above integral may be performed by 
hand or using symbolic manipulators. However, it is 
more convenient to use numerical integration in order 
to add flexibility to the programming. Nevertheless, 
some cases will not allow for hand or symbolic 
manipulators’ integration, thus it is common to find 
numerical integration routines associated with any 
finite element program. One of the popular techniques 
used for numerical integration is the Gauss-Fegendre 
quadrature which produces almost exact integration 
for polynomials that may reach up to 24 th order. 12 

2.3. 2-D Problems 

In 2-D problems, the interpolation polynomials should 
include combinations of the x and y coordinates. To 
identify the terms, the common practice is to select 
them from what is known as the Pascal triangle (See 
Figure 1). The number of terms selected should be 
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equal to the number of degrees of freedom and it is 
recommended that the terms should create a 
symmetric shape in that triangle. 


x y 

x A 2 xy y A 2 


x A 3 x A 2y 


xy*2 y A 3 


X A 4 x A 3y x A 2y A 2 x^3 y A 4 


Figure 1. Pascal triangle 

The Pascal triangle is quite handy in many problems, 
however, if you strict your work to using the full 
polynomials in both directions, then the generation of 
the terms becomes quite straight forward by 
multiplying both x and y-polynomials. For example, 
for a quadrilateral element with a single degree of 
freedom per node, we need a linear polynomial in 
each of the directions, hence, we may write: 

'I x' 

Lr xyj 

y)= [I x y xyJ 

Then we proceed with the H(x,y) vector just as in the 
1-D cases: 

\H(x u n) j' 

[H{x 2 ,y 2 )\ 


H (x,y)= 


-M 1 HI <y] 


[T]= 


From which we obtain: 

LA/(x,y)j=LH(x,y)J[r 1 ] 

Note, again, that the order of the polynomial terms in 
the H vector will not affect the resulting shape 
functions similar to what we had in the 1-D problems. 
At this point, we still have both problems when 
generating the higher order elements, namely, the 
matrix inversion and the numerical integration. In the 
following section, we will attempt generating the 
matrices using exact integration to avoid the 
numerical integration problem. 


3. EXACT INTEGRATION 

To avoid the errors introduced by the numerical 
integration, we may resolve to exact integration using 
symbolic manipulators (wxMaxima® , Mathematica® 
, Mable® , etc ...) but then we are back to the 
problem of rewriting the element matrices into our 
code, with all the problems that may include, or write 
the whole finite element model using the symbolic 
package which is normally extremely slow when it 
comes to numerical manipulations compare to 
numerical coding packages (Octave® , Matlab® , etc 
...) or programming languages (Fortran, C++, etc ...). 
Thus, we resolve to write an algorithm to perform the 
exact integration on the numerical package for 
specific elements with unknown number of degrees of 
freedom. This is where the H(x) row-matrix becomes 
very handy. Since the row-matrix is composed of 
simple polynomial terms, their differentiation, and 
later integration, is readily evaluated. 

Let’s examine the case of a bar element with n nodes 
(and n degrees of freedom). For the bar element, the 
element matrix, stiffness matrix, is evaluated using 
the integration: 

f 2 

k e =J EA{N’(x)}[N'(x)\dx 

For the sake of the illustration, we will assume that 
the modulus of elasticity and the cross-section area 
are constants, thus, we may write: 

Kr&f [rYfH'WHH'Wjir 1 ]* 

But, the transformation matrix is also constant, thus, it 
may be dragged out of the integration leaving us with 
the derivative of the H(x) vector. 

k= EAlT'fj' {hT(x)}{hT(x)\dxlT'] 

X , 

In this case, the T matrix may be given as: 

\H(x,)\' 

[7]_ [H{x 2 )\ 

[H(x n )\ 


If we examine the terms of the H(x) vector, we may write: 

[H(x)J= [l x x 2 ...J or /?,— x 1-1 /= 1,2,3,...,/? 

Thus the first derivative may be written as: 

[H’(x)\- [0 1 2x ...J 
or 

h’= 0 /=1 

' (/- 1)Y 2 /= 2,3 
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Hence, we may write the matrix: 


where 


{H'(x)}[H'(x)\=[G(x)} 


n-h'h'- 0 when /V j= 1 

9lJ ' J (i~ 1)(y - 1)^'" 2)+(r2) i,j= 2,3,...,n 

From that, we may get the integration: 

k= EA[T'] t \ 2 [G(x)]dx[T']= EA[TY[G][T'] 

X, 

where 

0 when /V j= 1 

9a= (i~ i)(y~ 1) / j/-2)+(;-2)+i_ j/-2)+(/-2)+i, i,j=2,3,...,n 

(/-2)+(y-2)+r 2 1 ' 

Similarly, for any 1-D problem involving any derivative, we may readily obtain the general term for the G 
matrix and perform the integration to obtain the G* matrix. 

Let’s examine the case of a beam element with n nodes (and 2n degrees of freedom). For the beam element, the 
element matrix, stiffness matrix, is evaluated using the integration: 




For the sake of the illustration, we will assume that the modulus of elasticity and the cross-section second 
moments of area are constants, thus, we may write: 


K= Elf l r’] T { H "(x))l H "(x)J[ r ’]d* 


But, the transformation matrix is also constant, thus, it may be dragged out of the integration leaving us with 
the derivative of the H(x) vector. 


K= E/[r 1 ] r j'{H' , (x)}[H"(x)Jdx[r l ] 


In this case, the T matrix may be given as: 


[T]= 


[H'W\ 

[H(x 2 )\ 

[H'(x 2 )\ 




[H(x n )\ 

H'(x n )\ 


If we examine the terms of the H(x) vector, we may write: 

[H(x)\= [l x x 2 ...J or h,= Y 1 /= 1,2,3 ,...,n 

Thus the second derivative may be written as: 

[H"(x)]= [0 0 2 6x ...J 
or 

h"= 0 /=1 ’ 2 

' (/- 1)(/- 2) x' -3 /= 3,4,.n 
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Thus, we may write the matrix: 

{H"(x)}[H"(x)\=[G(x)] 

where 

h " h " = 0 when /v J= 1 > 2 

yii ' ' (/— 1)(/- 2)(y— 1)( j- 2)x* /_3)+(7 “ 3) i ,j= 3,4,. ..,n 

From that, we may get the integration: 

k e = EI[T yf [G(x)]dx[r’]= E^irVlG'nr 1 ] 


where 


0 when /v j= 1,2 


9a= (/- 1)(/'-2)(y-1)(y-2), 


(/- 3)+(y- 3)+1 


(4' 


3)+(y- 3)+1_ ^/-3) + (y-3)+1^ 


i,j= 3,4,..., a? 


For 2-D problems, the problem becomes a little more complex. The general term of the H(x,y) vector will 
depend on the order at which we select the terms from the Pascal’s triangle. For illustration, let’s assume that 
we are going to generate a rectangular element for plate bending that ensures continuity of slope ( C 1 element). 
In this case, the number of nodes in the element may be selected to be n 2 where n is the number of nodes per 
side of the element (n=2,3,4, ...) in this case the number of degrees of freedom per element will be 4n 2 , and the 
x and y polynomials, each, will be of 2n-l order. We may write: 



1 1 

X 

A 2 . 

. x 2 ^ 1 

[H'(x,y)]= 

y 

xy 

A 2 / . 

x 2 ^ 1 y 


jT' 


x 2 /^ 1 ’ 

. x 2 ^ 1 /" 


Thus, we may write the H(x,y) vector as: 

[H(x, y)\= [1 x 


A 2 


rr 1 


y xy ... x 2 ^ 1 )/ 2 ^ 1 ] 


Whose general term may be given by: 

h k = Hjj= x j ~ 1 y- 1 where k= j+2n(i- 1), i,j= 1,2,. ..,2n 

And the derivatives of the H vector, may be written as: 

hk,xx= (j~ 1)(/-2)x y ~V _1 

hk,yy=(i- 1)('-2)x y_ V" 3 

hk,xy= (y -1 )( '“ i) x‘~ 2 y _ 2 

k= j+2n(i- 1) 

The element stiffness matrix for the static loading problem of thin plates may be written as: 

WT [TyiCfmcfT^dydx 


Y i 
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Where: 

I J | 

[H y {x v y,)\ 

[H xy {x u y,)\ 

[H{x 2 ,y 2 )\ 

[H x {x 2 ,y 2 )\ 

[H y {x 2 ,y 2 )\ 

[H Xy {x 2 , y 2 )\ 

[H x (x rf ,y rf }\ 

iHyiXrf, Yrf )\ 

^Jx^yJ J 


m 


4 rfx An 




And: 


[C] 


3x 4 rf~ 


'[H xx {x,y)\ 

[Hyy{x,y)\ 

2[H xy (x,y)\ 


Whose general term may be written as: 

^1 k~ ^k,xx 
C 2k = hk,yy 
°3k = K,xy 

k= 1,2,. ..,4/f 

[QJ is the plate stiffness which may be written, for the sake of compactness, in the form: 

Ta b O' 

[Q]= b a 0 
0 0c 


The matrix [G] may be written as: 


[G]= [C] r [Q][C] 


where the general term may evaluated by: 

Qkl~ XX^I XX+ yyhj y^+b(h/ < yyh/ XX+ xxhl yy^+A- Ch^ ^yhyxy 

Let’s now examine each term. 

Ur 1 )Ur 2)(y-1 )Ur 2)x^ 6 y^" 2 

h k ,yyh liX y={i k - 1)(/*-2)(y- 1)(y-2)^ 4 y^ 4 

h k ,„h liyjr ( j~ Wr 2 )(/,- 1 )(//- 2 )x^'- 4 y w '“ 4 

h k ,yyh Lyy =(i k - 1)(4-2)(/- 1)(/ r 2)x h+jr2 y k+ir6 
hk.xyh,,^ Ur i)(i/f-1) (yr i)(//-1) x y ‘ +yr 4 y iff,r 4 
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Where 


f 


j k = mod 


k -1 


\ 


h = 


\ 2 n +1 j 

k ~ h 


+ 1 


2n +1 


+ 1 


r 


j, = mod 


h = 


/-I 

2n +1 

l~h 


\ 


+ 1 


2 n +1 


+ 1 


Performing the integration, we get: 

xxdydx= 


cr,=f j h> „h „dvdx= Mk Oik 1M 2) ( X ' >+;, ~ 5 - x!‘ +;r5 “ 


*i 7i 


Uk + jl 5) (//;+/: 1) 




nyr- y,**' r1 ) 


Similarly: 


o' 2 =/ / h k yy h l xx dydx= 

x, y i 

a 2=j ( Kyyh l ,xxdydx= 


dk "0(4 2) ( /, 1) ( 7/ 2) ( y'c+y/- 3_ X A+ yr 3y y>/3_ y> i,- 3^ 


x, y, ( jk+jl 3 ){ik+h 3) 


(4 "0(4 2) (7/ 1) ( 7/ 2) / „j k +jr 3_ v ,y ( ,+ 7-3^y>/-3_ y t +/;-3^ 


-{4 


*i yi 


h k xv hi Xv dydx= 

x, y, ’ {jk+jl 3)(//(+// 3) 


( jk+jl 3) (4+// 3) 

— (*2 

( jk+jl 1)(//(+// 5) 

1 '(V 


_ f f h h ^4 ^) (4 3) (// 1)(// 3) i_ 1 w 5_ ./>/- 5\ 

“4 , (/„+/,-1)(/;+/,-5) ( 2 1 My2 y ' 

a.rr' 


\ j k / ( 4 "^ ) (■// "^ ) ( 4 / / \/jk+ ji 3 v jk+ji 3y y*+ i 3 J k +ii 3^ 


y r ,-3) 


From which, we may get the integral of the general term of the [G] matrix in the form: 


g u dydx= g M = a{a ,+a 4 )+b{a 2 +a 3 )+4 Cor, 


ff 

*1 Xi 

In both the 1-D and 2-D problems illustrated above, we were able to create the element matrices exactly using 
the trick of integrating the 77 vector and its derivatives. However, the limiting condition on generating high 
order elements will always be the inversion of the T matrix which will become ill-conditioned with higher order 
elements. For example, the plate bending element described above worked only until the number of nodes 
became 16 (64 DOF) (7 th order polynomial in each of the x and y directions) using double precision numbers 
on the Octave numerical manipulator. Hence, a search for another technique was required. 

4. LAGRANGE POLYNOMIALS 
4.1. 1-D problems 

The Lagrange interpolation polynomial for an (n-l) th order polynomial, which may be used when n degrees of 
freedom are give, may be written in the form: 

|M 

4 


n n -1 


nx)-jfn^-mx) n 2 ( X ) ... w„(x)j 

i^l /= 1 X, X; 


pi 


f. 
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Where: 


n- 1 


X- Xj 


w,. ( x)= n — 


>1 *r x , 

j*i 


Thus, we may obtain the first derivative in the form: 

n- 1 

N';(X) = 


n~1 H n~ 1 y— y 

'i(x)=j —— n — l 

fei x r xj 

kt i k 


i*k 

i*i 

We may also obtain further derivatives for use in the element equation, however, let’s remember that the 
function interpolation is based on the values of the function, thus, if we need to have any of the derivatives of 
the function as a degree of freedom, we cannot use this interpolation method. Nevertheless, we have obtained 
the shape functions without having to go through the matrix inversion problem described in the previous 
section. On the other hand, performing the exact integration will not be readily available for the general term 
since it now involves multiplication of several linear terms. Luckily, as mentioned earlier, the Gauss quadrature 
integration provides high degree of accuracy for a relatively high order polynomial integration. 

Note that he order of the Lagrange polynomials used above has to coincide with the order of the node 
numbering unlike the cases we described in sections 2 and 3 above. 


4.2. 2-D Problems 

With the same advantages and disadvantages of the 1-D problem, the derivation of the shape functions for a 2- 
D problem comes straight forward by multiplying the polynomials in the x and y directions. 


f ^^lP ! U^t\jFvr iN " {x ’ y) w,2<x ' y) - w « (x ' r)l 


f 12 


f „ 




In a vector, we usually use a single index to point to the different element, thus, we need to create a relation 
between the vector index and the i and j counters of the nodes counters in a 2-D element. A simple relation for 
counting the nodes in the x-direction first would be in the form: 

m=i+{j- 1 )n x 
i= 1,2 ,...,n x 
7= 1,2 ,...,n y 
m= 1,2,. ..,n x x n y 

Thus we may write the general shape function in the form: 


yr y> 


k£ i 


l*j 


For which we may evaluate the first derivatives to be: 


n -1 


N n 


n - 1 x- x u y- y, 


<*.y)=/Mx,y)=X—n^Q- 


N, 


P=1 I P K = 1 
p? i k? p 

kt i 

Ry~ 1 yr_ y, 1 

A(r 


// 


!*i 


n„- 1 


x _ x £ y _y^ p y _ 

kt i q* j lit q 


y- y_ 

Yi 


>*q 

i*i 


n- 1 


Pr 1 X~ X ^ 


n v - 1 


N m Xy (x,y)= Ny ^{x, y)= T-fl-- Y-fl — 

x-x p U x-x k fa y-y q U y- 

pt i k£ p qt j It q 

kt i 


y- y_ 

y 


l*q 

l*j 
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5. MODIFYING THE INTERPOLATION FUNCTION 

Borrowing from the Spline interpolation techniques, we may work around with the Lagrange polynomials to 
obtain polynomials that may include derivatives of the function. In the following, we will work with functions 
that are required to satisfy the value of the function and its first derivative at each node of the element. We will 
start be deriving the cubic polynomial that will fit 2 values and two slopes, then, we will describe the process 
for quintic polynomial, followed by the general approach for obtaining the (2n-l) th polynomial that may be 
used for n-node elements. 


5.1. The Polynomials 

Using the similarities we obtained in the 3 rd and 5 th order polynomial derivation, we will present a 
generalization in the following section. We may write the general relation for a function of the (2n-l) th order 
that uses n pairs of polynomials as: 


f(x)=y sax) 


Where: 


S(x)= (C,+ C 2 (x- x,))f| (j-pit 
This function needs to satisfy the conditions: 


i*i 


$(*/)= f, 

s',(x,)=r, 

For the first condition, we get: 

C 1= f i 

To proceed with the second condition, we get the slopes as: 


s:w=c 2 p(^) + 2(c, + c 

Mi 


n 


k± i 


X- x k 


x k ) '= 


p (^) 2 


i*i 

i*k 


Substituting: 


Which gives: 


s',(x,)=r,= c 2+ 2c,| -J- 


k£i 


C 2 = f 2 f, 


fa x r x k 

k+i 


Thus we may write the pair of polynomials as: 

Si(x) 


f, 


fa x r x k 

kt i 


V ")0 <s)' 


Rearranging, we obtain: 


S(x)= ^) + r ,(x-x;.,p 

\ kti jn pi 


From which we may write the shape functions as: 

Ni,i{x)= /l- 2(x 


N l 2 (x)= {x- x ( ) 



pi 

i= 1,2 
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In 2-D problems, the interpolation functions will be a straight forward generalization of the above procedure 
giving the function in the form: 

f(x,y)=f_ S m (x,y) 


m= 1 

n x x ny 


Where: 


Smi x )~ (£ 1 +C 2 (x x,)+C 3 (y y y )+C 4 (x X/)(y 


«'CI (W 7 ,j 

kV i i ? j 


m= i+{j~ 1)* n x 
i= 1,2 ,...,n x 
j= 1,2 ,...,n y 
m= 1,2,. ..,n 

This function needs to satisfy the conditions: 

^m[Xj, yj)= f ij 
$n,x{Xi,yj)= fij.x 

S m , y (*,,y,)= fij,y 

^m,xy(Xj, yj)— f ij xy 

Which will result in a quartet of shape functions in the form: 

r v 

1 1 


N m ,x( x >y)= 


\-2(x- x, )Y, 


*=i x i - x k 
k^i J 


P =i yj - y P 

p*j j 


X 

n 

/=i 
l^i 


x — x, 


V X i ~ X l J 


j 

n 

9=1 

9*j 


y-y 9 

~y<i j 


n , \ n, , .2 n v , .2 

\ kti / I* i qtj 


N rr Jx, y)= (X- x ( ) |l-2(y- y y ) £ -~ 
N m , 3 (x,y)=^ 1- 


/#/ 

/= 1,2, ...,n x 
j= 1,2 ,...,n y 
m= 1,2,...,n x x n y 

5.2. The Derivatives 

In the problems that make use of the shape functions described in the previous section, the finite element model 
usually requires the first and second derivatives of the shape function to perform the required calculations. In 
this section, we will present those derivatives for the 1-D problem without elaborations. 


N/,,(x)= /1-2(x-x i )| 

\ ^ 

N i 2 (x)=(x- x,] 



j*i 


A/'„i(x)=/- 


l fei X- X k f j= \\x r XjJ I 

\ kti / j* / \ 


1 - 2(x- x,V 


X- X, 


to* i 


\fo x ~ x ' n 

/x-x,\ 2 

Fi "(x r x ; ) 2 U 

V x -x ; 7 


f Vi 


i*i 

i*i 
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6. NUMERICAL RESULTS 

All the test cases used were confined to generating the element stiffness matrix for different problems for 
isoparametric C° and C 1 continuous elements. The number of nodes were increased and the Eigenvalues of the 
matrix were evaluated until the lowest (negative) Eigenvalue reached an absolute values larger than 10 4 at 
which the element was considered a failure since the matrices investigated were supposed to be positive semi- 
definite. The results for the different problems are presented below. 

Table 1 presents the change of the values of the minimum Eigenvalue of the stiffness matrix generated for a bar 
element. If we ignore the numbers that have an absolute value less than 10 10 considering them as a numerical 
zero, we may still find that increasing the number of nodes per element start to introduce negative Eigenvalues 
until we reach the failure criterion set in this work as 10‘ 4 . Note that, unpredictably, the exact integration 
method failed before the numerical integration at 15 nodes per element. On the other hand, the Lagrange 
polynomial method, which does not use the transformation matrix, was able to produce satisfactory results up 
to 23 nodes per element which approaches the limit of the 12 point numerical integration technique used. 

Table 2 presents the results obtained for a beam element. Again we find that the exact integration and numerical 
integration elements, both, failed around the 18 and 16 DOF respectively which coincides with the number of 
DOF at which failure occurred in the bar element. Meanwhile, the Modified Lagrange method was able to 
continue all the way up to 24 DOF which reaches the limit imposed by the numerical integrator. 


Table 1. Change of the lowest Eigenvalue for the Bar problem using different methods 


Number of Nodes 

Classical Method 

Exact Integration 

Lagrange Polynomials 

2 

0.0000 

0.0000 

0.0000 

5 

0.0000 

1.3550 * 10- 15 

-3.2092 * 10 14 

13 

-4.3771 * 10' 11 

1.8269 * 10' 8 

-4.7237 * 10" 9 

14 

-3.5311 * 10' 10 

1.2446 * 10' 6 

-8.3456 * 10' 9 

15 

-6.7339 * 10' 9 

-7.2761 * 10 3 - Failure 

-1.7247 * 10' 8 

16 

-4.0141 * 10" 5 


-6.8237 * 10~ 8 

17 

-0.0063877 - Failure 


-2.7224 * 10‘ 7 

24 



-0.000284-Almost! 
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Table 2. Change of the lowest Eigenvalue for the Beam problem using different methods 


Number of Nodes (DOF) 

Classical Method 

Exact Integration 

Modified Lagrange Polynomials 

2(4) 

.1.8447 * 10' 15 

-1.6821 * 10' 15 

-4.0007 * 10 16 

5(10) 

-2.7423 * 10' 10 

-2.3674 * 10' 10 

-7.2653 * 10~ 13 

6(12) 

-2.7096 * 10~ 8 

-4.4350 * lO' 10 

-3.3600 * 10 10 

7(14) 

9.5001 * 10" 8 

4.6870 * lO' 8 

-1.4082 * 10' 9 

8(16) 

-2,4339 - Failed 

-2.7862 * 10~ 4 

-1,0945 * 10' 9 

9(18) 


-0.045658 - Failed 

-3.2284 * 10' 8 

11(22) 



-7.7276 * 10' 7 

12 (24) 



-2.9915 * 10' 5 


Table 3. Change of the lowest Eigenvalue for the Plate problem using different methods 


Number of Nodes 
(DOF) 

Classical Method 

Exact Integration 

Modified Lagrange 
Polynomials 

4(16) 

-2.4664 * 10' 16 

-9.9632 * 10" 16 

-4.2013 * 10 15 

9(36) 

-4.0203 * 10' 13 

-4.3150 * 10" 13 

-1.1133 * 10 14 

16 (64) 

-2.1111 * 10' 9 

-3.6665 * 10 10 

9.6377 * 10' 16 

25(100) 

-2.6320 * 10 +5 - 
Failed 

-1.9562 * 10 +7 - 
Failed 

-2.7317 * 10 14 

100(400) 



-4.5663 * lO’ 6 

121(484) 



-2.1076 * 10‘ 4 -Almost! 

144 (576) 



-0.023264 - Failed 


Table 3 presents the results for a plate bending problem. The same patterns observed in the bar and beam 
elements can be observed in the plate problem for the modified Lagrange case. 


7. CONCLUSIONS 

Several problems were examined in this work to 
demonstrate the limitations imposed on the creation of 
super elements with high number of nodes and 
degrees of freedom. The element matrices were 
generated using simple, Lagrange, and modified 
Lagrange polynomials and the integration was 
performed numerically using Gauss-Lagrange 
quadrature or exact when possible. Problems 
involving 1 -D elements with 1 or 2 DOF per node and 
2-D problems with 4 DOF per node were examined 
by increasing the number of nodes per element until 
the lowest (negative) Eigen value of the stiffness 
matrix had an absolute value greater than 10 4 . From 
the results obtained above we may conclude the 
following: 

> Automated generation of higher order 
isoperimetric C° and C 1 continuous element model 
was enabled using an open source numerical 
manipulator (Octave®) using the different 
methods presented in this work. 

> Generating higher order elements is 
computationally expensive compared to lower 
ones. Thus, the h-version is more efficient for the 
same accuracy than the p-version. 


> Using numerical integration to evaluate the 
element matrices is efficient and accurate enough 
for most practical purposes 

> Using Lagrange and modified Lagrange 

polynomial increases the stability of the results for 
higher order element problems when they become 
a necessity to use. In such cases, the limitation on 
the results’ accuracy will be imposed by the 
accuracy of the numerical integration scheme. 

Further work with the modified Lagrange 
polynomials, proposed in this work, is needed to 
enable the inclusion of higher order derivatives in the 
interpolation polynomials. Such research will help in 
the generation of elements with the needed derivatives 
without falling into the trap of the ill-conditioned 
transformation matrix. 
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